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Dendrites are the primary component of solidification microstructures in metals. Their properties 
have been a topic of intense study in the past 10-15 years. Experiments [1,2] by Glicksman and 
coworkers on succinonitrile (SCN) and other transparent analogues of metals have provided tests 
of theories of dendritic growth, and have stimulated considerable theoretical progress [3-5]. The 
experiments have clearly demonstrated that naturally growing dendrites possess a unique steady 
state tip, characterized by its velocity, radius of curvature and shape, which leads the to a time- 
dependent sidebranched dendrite as it propagates. 

Insight into the steady state dendrite problem was first obtained from local models [6-9] describing 
the evolution of the interface, and incorporating the features of the bulk phases into the governing 
equation of motion for the interface. These models showed that a nonzero dendrite velocity is 
obtained only if a source of anisotropy - for example, anisotropic interfacial energy - is present in 
the description of dendritic evolution. It was then shown that the spectrum of allowed steady state 
velocities is discrete, not continuous, and the role of anisotropy was understood theoretically, both 
in the local models and the full moving boundary problem [5,10,11]. Moreover, only the fastest 
of a spectrum of steady state velocities is stable, thus forming the operating state of the dendrite. 
It is widely believed that sidebranching is generated by thermal or other statistical fluctuations on 
a microscopic scale, which are amplified by advective diffusion. This body of theoretical work is 
generally known as solvability theory. 

Brute force solution of the time-dependent Stefan problem requires front tracking and lattice de- 
formation to contain the interface at predefined locations on the grid [12]. The phase-field model 
avoids this problem by introducing an auxiliary continuous order parameter <f>( r) that couples to 
the evolution of the thermal field. The phase field interpolates between the solid and liquid phases, 
attaining two different constant values in either phase, with a rapid transition region in the vicinity 
of the solidification front. The level set of <f>( r) = 0 is identified with the solidification front, and 
the dynamics of (j) are designed to follow the evolving solidification front [13-23]. The phase-field 
parameters can be derived from parameters of the Stefan problem [14,24], however this mapping 
is not very sensitive to the precise form of the phase-field model [25]. 

While the phase-field model finesses the problem of front tracking, it is still prohibitively expensive 
for large systems, because the grid spacing must be small enough everywhere that the phase-field 
model converges to the the sharp interface limit [14,24]. Caginalp and Chen [26] showed rigorously 
that the phase field model converges to the sharp interface limit when the interface width (and hence 
the grid spacing) is much smaller than the capillary length. This result is necessary for acceptance 
of the phase field model, but is not sufficient for computational tractability in the experimentally 
relevant regime. 

However, more recently Karma and Rappel [24] presented a different asymptotic analysis in powers 
of the ratio of the interface width to the diffusion length. Their procedure allows selection of 
parameters such that the phase field model corresponds to the sharp interface limit when the 
interface width (and hence the grid spacing) is of order the capillary length - a much more tractable 
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regime. Furthermore, their improved analysis allows the kinetic coefficient to be tuned to zero, 
which corresponds to the experimentally realized situation at low undercooling in succinonitrile 
[2], Karma and Rappel’s numerical results are in excellent agreement with solvability theory at 
dimensionless undercoolings as low as 0.30, but fail to access the range of experimentally- realizable 
undercoolings near 0.1. What is needed is an effective adaptive technique [27] which dynamically 
coarsens the grid spacing away from the front. 

In this work we show how the phase field model can be solved in a computationally efficient manner 
that opens a new large-scale simulational window on solidification physics. Our method uses a finite 
element, adaptive-grid formulation, and exploits the fact that the phase and temperature fields 
vary significantly only near the interface. We illustrate how our method allows efficient simulation 
of phase-field models in very large systems, and verify the predictions of solvability theory at 
intermediate undercooling. We then present new results at low undercoolings that suggest that 
solvability theory may not give the correct tip speed in that regime. 

We model solidification using the phase-field model used by Karma and Rappel [24]. We rescale 
temperature T by U = cp(T — Tm)/L , where cp is the specific heat at constant pressure, L is the 
latent heat of fusion and Tu is the melting temperature. The order parameter is defined by <£, 
with 0 = 1 in the solid, and </> = — 1 in the liquid. The interface is defined by (ft = 0. We rescale 
time by r G , a time characterizing atomic movement in the interface, and length by a length 
characterizing the liquid-solid interface. The model is given by 


A 2 (fi) 
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where D = ar 0 /W% and a is the thermal diffusivity, and where A controls the coupling of U 
and 0. Anisotropy has been introduced in Eqs. (1) by defining the width of the interface to be 
W(n) = W 0 A(n) and the characteristic time by r(n) = r 0 A 2 (fi) [24], with A(n) £ [0,1], and 


A(n) = (l-3e) [l + A 


The vector n = (<t>, x x + <t>, y y)/(<f> 2 x + <t> 2 y ) 1/2 is the normal 
to the contours of </>, and and c p }V represent partial derivatives with respect to x and y. The 
constant e parameterizes the deviation of W (n) from W Q - We expect the results to be similar for 
other definitions of anisotropy [14]. 


We use the asymptotic relationships given in [24] to select the parameters in Eqs. (1) such that it 
operates in the sharp interface limit, defined by U at the interface satisfying U\ n t = —d(n)K— j3(n)V n . 
The variable d(n) is the capillary length, k is the local curvature, /?(n) is the interface attachment 
kinetic coefficient and V n the normal speed of the interface, all in dimensionless form. In terms 
of A(n), d{n) = d 0 [A(n) + c^A(n)], where d 0 = 0.8839/A and 6 is the angle between n and the 
x-axis. In this formulation, the constants W 0 , r 0 and A may be chosen so as to simulate arbitrary 
values of (3 . In particular, A = 1.5957£> makes (3 = 0 [24], a limit which is appropriate for SCN. 


We compute four-fold symmetric dendrites in a quarter-infinite space, initiated by a small quarter 
disk of radius R 0 centered at the origin. The order parameter is initially set to its equilibrium value 
— — tanh((|r| — R 0 )/y/ 2) along the interface. The initial temperature decays exponentially 
from U — 0 at the interface to — A as x — > oo. 


We simulate Eqs. (1) on an adaptive grid of linear isoparametric quadrilateral and triangular finite 
elements, formulated using Galerkin’s method. Following Ref. [28], elements are arranged on a 
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two dimensional quadtree data structure, which makes our code scalable when implemented using 
dynamic memory allocation. The largest system sizes we have considered thus far correspond to 
2D uniform lattices having 2 17 x 2 17 grid points. The grid is locally refined to have a higher density 
of elements in the vicinity of the interface, identified by large fluxes in a composite field based on 
both <p and U. Typically, the grid is adapted about every 100 time steps, which permits <f> and U 
to remain within the refined range between regriddings. We allow a difference of at most one level 
of refinement between neighboring quadrilateral elements. In such a case the quadrilateral element 
of lower level of refinement has an extra side node. The extra nodes are resolved with triangular 
elements. 

On an adaptive grid, the concept of a grid spacing is replaced by that of a minimum grid spacing 
Ax m in, representing the finest level of spatial resolution. We found that for solutions to converge 
properly, the grid must be layered such that the highest density of elements appears around the <p 
interface, while the [/-field, whose width is of order D/V n can be encompassed by a mesh of larger 
grid spacing. Convergence of our solutions is relatively insensitive to A:r m j n . For a test case of 
dendrites grown at A = 0.55, D — 2, c = 0.05, and integration time step dt = 0.016, our solutions 
for the steady state velocity converge to that given by solvability theory to within a few percent 
for 0.3 < Ai m j n < 1.6. 

Fig. 1 shows a dendrite 10 5 time steps into its evolution computed using our adaptive grid method, 
using the parameters just mentioned above. The system size is 800 x 800, with Ax m j n = 0.78, 
and about half of the computational domain is shown. Sidebranching is evident, and arises due 
to numerical noise. This calculation took approximately 10 cpu-hours on a Sun UltraSPARC 2200 
workstation. 

We examined the cpu-scalability of our algorithm with system size by growing dendrites in systems 
of various linear dimension L B and measuring the cpu time Rf for the dendrite branches to traverse 
the entire system. We once again use the same parameters as above, except Ax m j n = 0.4. The 
relationship between Rf and Lg is shown in Fig. 2, where we see that Rf ~ Tg. The number of 
calculations performed, per time step, is proportional to the number of elements in our grid, which 
is set by the arclength of the interface being simulated multiplied by the diffusion length D/V n . For 
a parabolic shape the arclength ~ Lg. Thus, since the dendrite tip moves at a constant velocity V n , 
then Rf = [RfDV 2 / &x 2 m ]L 2 B , where Rf is a constant that depends on the implementation. The 
cpu time Rf needed to compute the same case on a uniform grid scales as Rf = [/[“/(UnAx^jjL^. 
For large system sizes, Rf / R^ ~ L B . 

We tested the effective anisotropy of our dynamically adapting lattice in two ways. Following the 
method outlined by Karma [24], we find an equilibrium shape for the interface when the background 
field is adjusted dynamically so as to maintain the velocity of the interface at zero. The effective 
anisotropy is inferred by fitting an equation to the computed interface. We found c e ff to be within 
5% of the intended value for input e = 0.02 - 0.04. We also tested for grid anisotropy by rotating 
the grid by 45 degrees, which should represent the lowest accuracy for square elements. In this 
case, the steady state tip-velocity was within 1% of its value in the original orientation. 

We further verified our algorithm by comparing measured tip-velocities and shapes for dendrites 
grown using the same undercoolings, parameter sets and systems sizes reported in [24], We found 
very good agreement for A = 0.65,0.55,0.45,0.30. ^Ve next investigated the effect of system size. 
Fig. 3 shows the time evolution of tip-velocity for several undercoolings and system dimensions. 
The two cases for A = 0.65 are typical of results at intermediate A, showing a relatively rapid 
leveling to an asymptotic speed within a few percent of that predicted by solvability theory. 

At lower A, however, we found that the tip-velocity deviates from that predicted by solvability 
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theory. Fig. 3 also shows the evolution of the tip-velocity for A = 0.25 in two different sized boxes. 
Whereas the computed tip- velocity falls a few percent below the solvability value in the 6400 x 400 
box, it exceeds by 8% the solvability value in the 6400 x 3200 box. This effect is even larger at 
A = 0.1, also shown in Fig. 3, where the tip speed is about 3 times larger than that predicted by 
solvability theory. 

The explanation for this behavior is that at low A, the thermal fields of the two dendrite branches 
overlap, violating the assumptions of solvability theory, which models an isolated single dendrite. 
At large undercooling, each dendrite arm quickly outruns the other’s thermal boundary layer, and 
solvability theory should apply. (See Fig. 1, A = 0.65.) The conditions of solvability theory can 
also be approximated at lower undercooling if simulations are performed in a domain which is 
small in one direction. For the simulation performed with A = 0.25 in the small box (6400 x 400), 
the branch in the y-direction is extinguished by its interaction with the wall and agreement with 
solvability theory is obtained. However, when both branches are present, as in the simulation with 
A = 0.25 in the larger box (6400 x 3200), their interaction leads to an increased tip- velocity because 
the dendrites are embedded in a circular rather than parabolic diffusion field. This is clearly seen 
in Fig. 4, where the dendrite shape and its associated field are shown for A = 0.10 ( D = 13, 
do = 0.043, e = 0.05, Ax = 0.78, dt = 0.08). The dendrite arms never became free of each other in 
this simulation, causing the observed deviation from solvability theory shown in Fig. 3. This latter 
simulation was performed in a 102400 x 51200 domain, chosen to contain about 10 D/V n . We note 
that the ratio of the largest to smallest element size in this simulation is 2 17 . A fixed mesh having 
the same resolution would contain 9 x 10 9 grid points, clearly beyond current computing capability. 

We can estimate the time t* when the growth of the dendrite tip crosses over from the transient 
regime where the branches interact to that where they become independent by equating the length 
of the full diffusion field, 3(Z?t*)'^, to the length of a dendrite arm, V n t*. This gives the crossover 
time as t* = 9D/V . „ 2 . The values for t* corresponding to the cases A = 0.65, 0.25 and 0.10 in Fig. 3 
are 2.5 x 10 3 , 1.6 x 10 4 and 5.9 x 10 7 , respectively. Inspection of Fig. 3 confirms this scaling. 

These results have important implications when comparing theory to experimental observations at 
low undercooling. We find that in this regime, the appropriate theory to use is one which explicitly 
takes into account the long range effects of other branches [29]. In particular, study of real dendrites 
with sidebranches, growing at low undercooling, will require such treatment. An investigation of 
this effect, as well as results in 3D, applications to directional solidification and other solidification 
processes, and a more detailed description of our algorithm will appear in future publications. 

We thank Wouter-Jan Rappel for providing the Green’s function steady-state code used to test some 
of our simulations, and Alain Karma for generously providing us with his unpublished results. 
We also thank Robert Almgren and Alain Karma for helpful discussions of our results at low 
undercooling. This work has been supported by the NASA Microgravity Research Program, under 
Grant NAG8-1249. 
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Figure 1: A dendrite grown using the adaptive- 
grid method for A = 0.55, D = 2, e = 0.05. 
Clockwise, beginning at the upper right the fig- 
ures show contours of the {/-field, the contour 
<j> = 0, contours of the <^-field and the current 
mesh. 



Figure 3: The time evolution of the tip-velocity 
for undercooling A = 0.65,0.25 and 0.10. 
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Figure 2: CPU time vs. the system size, illustrat- 
ing the computing time for a dendrite to move 
through the system of linear dimension L n using 
our adaptive mesh method. 



Figure 4: Dendrite silhouette and isotherms from 
-0.01 to -0.9 for undercooling AT = 0.1. Full 
domain dimensions are 102,400 x 51,200. The 
dendrite tip is approximately 1,300 units from 
the origin, while the temperature field spreads 
to about 5, 000 units. 
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